function outputdata = make3Dgraphs(paramA,paramB,paramC,AZIM,ELEV)

% Set universal view point for 3D graphs. Can be overridden via input.
global AZ;
global EL;
if nargin<5
    AZ = -37.5;
    EL = 30;
else
    AZ = AZIM;
    EL = ELEV;
end


clc


% start export output to txt
% delete([param1 '_' param2 '_' param3 '.txt']); % disabled automatic deletion of old
% logs for other users to avoid unwanted side effects - HM

%diary(['Output/' paramA '_' paramB '_' paramC '.txt']); % disabled diary function. Enable if
% needed. See also diabled diary('off') line at the end of this file. - HM






% masks without OC
fullmask = 3;
graymask = 4;
whitemask = 5; % 2 = WM inckuding the optical chiasm

%% full slice
fig=figure(101);
set(fig,'Position',[0 0 1000 1000])
disp(['2D full fit for ' paramA ' and ' paramC]);
%fig bottom left
temp=zeros(10,10);
temp(1:4,8:10)=1;
subplotpos=find(temp==1);
[outputdata.fullslice_A_C.x,...
    outputdata.fullslice_A_C.y,...
    outputdata.fullslice_A_C.z,...
    outputdata.fullslice_A_C.regressionresult,...
    outputdata.fullslice_A_C.x_unprocessed,...
    outputdata.fullslice_A_C.y_unprocessed,...
    outputdata.fullslice_A_C.z_unprocessed,...
    outputdata.fullslice_A_C.GMmask,...
    outputdata.fullslice_A_C.WMmask] = marmplots(paramA,paramC,[],fullmask,101,subplotpos);

disp(['2D full fit for ' paramB ' and ' paramC]);
%fig bottom right
temp=zeros(10,10);
temp(7:10,8:10)=1;
subplotpos=find(temp==1);
[outputdata.fullslice_B_C.x,...
    outputdata.fullslice_B_C.y,...
    outputdata.fullslice_B_C.z,...
    outputdata.fullslice_B_C.regressionresult,...
    outputdata.fullslice_B_C.x_unprocessed,...
    outputdata.fullslice_B_C.y_unprocessed,...
    outputdata.fullslice_B_C.z_unprocessed,...
    outputdata.fullslice_B_C.GMmask,...
    outputdata.fullslice_B_C.WMmask] = marmplots(paramB,paramC,[],fullmask,101,subplotpos);

disp(['3D full fit for ' paramA ' and ' paramB ' and ' paramC]);
%fig top
temp=zeros(10,10);
temp(2:9,1:6)=1;
subplotpos=find(temp==1);
[outputdata.fullslice_A_B_C.x,...
    outputdata.fullslice_A_B_C.y,...
    outputdata.fullslice_A_B_C.z,...
    outputdata.fullslice_A_B_C.regressionresult,...
    outputdata.fullslice_A_B_C.x_unprocessed,...
    outputdata.fullslice_A_B_C.y_unprocessed,...
    outputdata.fullslice_A_B_C.z_unprocessed,...
    outputdata.fullslice_A_B_C.GMmask,...
    outputdata.fullslice_A_B_C.WMmask] = marmplots(paramA,paramB,paramC,fullmask,101,subplotpos);

set(findall(fig,'-property','FontSize'),'FontSize',12)

print(['Output/' paramA '_' paramB '_' paramC '_fullslice'],'-dpng','-r300');

Npoints = outputdata.fullslice_A_B_C.regressionresult.fres.NumObservations;
SSEx1 = outputdata.fullslice_A_C.regressionresult.fres.SSE;
Nparax1 = outputdata.fullslice_A_C.regressionresult.fres.NumCoeffs;
SSEx2 = outputdata.fullslice_B_C.regressionresult.fres.SSE;
Nparax2 = outputdata.fullslice_B_C.regressionresult.fres.NumCoeffs;
SSEbiv = outputdata.fullslice_A_B_C.regressionresult.fres.SSE;
Nparabiv = outputdata.fullslice_A_B_C.regressionresult.fres.NumCoeffs;

disp(' ');
outputdata.fullslice_A_B_C.regressionresult.fres
outputdata.fullslice_A_C.regressionresult.fres
outputdata.fullslice_B_C.regressionresult.fres
disp(' ');

disp('Testing univariate fit object fullslice_A_C (Fit 1) against bivariate fit object fullslice_A_B_C (Fit 2) with a significance threshold of 0.001:');
[Fbivx1,pbivx1]=ftest([SSEx1 SSEbiv],[Nparax1 Nparabiv],Npoints);
disp(' ');
disp('Testing univariate fit object fullslice_B_C (Fit 1) against bivariate fit object fullslice_A_B_C (Fit 2) with a significance threshold of 0.001:');
[Fbivx2,pbivx2]=ftest([SSEx2 SSEbiv],[Nparax2 Nparabiv],Npoints);


disp(' ');
disp('The F-test for the univariate fit object fullslice_A_C vs. the bivariate fit object fullslice_A_B_C had the inputs');
disp('Nparams_univariate; Nparams_bivariate; Nsamples; SSE_univariate; SSE_bivariate of:');
disp([num2str(Nparax1) '; ' num2str(Nparabiv) '; ' num2str(Npoints) '; ' num2str(SSEx1) '; ' num2str(SSEbiv)]);
disp(' ');
disp('The F-test for the univariate fit object fullslice_B_C vs. the bivariate fit object fullslice_A_B_C had the inputs');
disp('Nparams_univariate; Nparams_bivariate; Nsamples; SSE_univariate; SSE_bivariate of:');
disp([num2str(Nparax2) '; ' num2str(Nparabiv) '; ' num2str(Npoints) '; ' num2str(SSEx2) '; ' num2str(SSEbiv)]);
disp(' ');disp(' ');disp(' ');

%% GM
fig=figure(102);
set(fig,'Position',[0 0 1000 1000])

disp(['2D GM fit for ' paramA ' and ' paramC]);
%fig bottom left
temp=zeros(10,10);
temp(1:4,8:10)=1;
subplotpos=find(temp==1);
[outputdata.GM_A_C.x,...
    outputdata.GM_A_C.y,...
    outputdata.GM_A_C.z,...
    outputdata.GM_A_C.regressionresult,...
    outputdata.GM_A_C.x_unprocessed,...
    outputdata.GM_A_C.y_unprocessed,...
    outputdata.GM_A_C.z_unprocessed,...
    outputdata.GM_A_C.GMmask,...
    outputdata.GM_A_C.WMmask] = marmplots(paramA,paramC,[],graymask,102,subplotpos);

disp(['2D GM fit for ' paramB ' and ' paramC]);
%fig bottom right
temp=zeros(10,10);
temp(7:10,8:10)=1;
subplotpos=find(temp==1);
[outputdata.GM_B_C.x,...
    outputdata.GM_B_C.y,...
    outputdata.GM_B_C.z,...
    outputdata.GM_B_C.regressionresult,...
    outputdata.GM_B_C.x_unprocessed,...
    outputdata.GM_B_C.y_unprocessed,...
    outputdata.GM_B_C.z_unprocessed,...
    outputdata.GM_B_C.GMmask,...
    outputdata.GM_B_C.WMmask] = marmplots(paramB,paramC,[],graymask,102,subplotpos);

disp(['3D GM fit for ' paramA ' and ' paramB ' and ' paramC]);
%fig top
temp=zeros(10,10);
temp(2:9,1:6)=1;
subplotpos=find(temp==1);
[outputdata.GM_A_B_C.x,...
    outputdata.GM_A_B_C.y,...
    outputdata.GM_A_B_C.z,...
    outputdata.GM_A_B_C.regressionresult,...
    outputdata.GM_A_B_C.x_unprocessed,...
    outputdata.GM_A_B_C.y_unprocessed,...
    outputdata.GM_A_B_C.z_unprocessed,...
    outputdata.GM_A_B_C.GMmask,...
    outputdata.GM_A_B_C.WMmask] = marmplots(paramA,paramB,paramC,graymask,102,subplotpos);

set(findall(fig,'-property','FontSize'),'FontSize',12)

print(['Output/' paramA '_' paramB '_' paramC '_GM'],'-dpng','-r300');

Npoints = outputdata.GM_A_B_C.regressionresult.fres.NumObservations;
SSEx1 = outputdata.GM_A_C.regressionresult.fres.SSE;
Nparax1 = outputdata.GM_A_C.regressionresult.fres.NumCoeffs;
SSEx2 = outputdata.GM_B_C.regressionresult.fres.SSE;
Nparax2 = outputdata.GM_B_C.regressionresult.fres.NumCoeffs;
SSEbiv = outputdata.GM_A_B_C.regressionresult.fres.SSE;
Nparabiv = outputdata.GM_A_B_C.regressionresult.fres.NumCoeffs;

disp(' ');
outputdata.GM_A_B_C.regressionresult.fres
outputdata.GM_A_C.regressionresult.fres
outputdata.GM_B_C.regressionresult.fres
disp(' ');

disp('Testing univariate fit object GM_A_C (Fit 1) against bivariate fit object GM_A_B_C (Fit 2) with a significance threshold of 0.001:');
[Fbivx1,pbivx1]=ftest([SSEx1 SSEbiv],[Nparax1 Nparabiv],Npoints);
disp(' ');
disp('Testing univariate fit object GM_B_C (Fit 1) against bivariate fit object GM_A_B_C (Fit 2) with a significance threshold of 0.001:');
[Fbivx2,pbivx2]=ftest([SSEx2 SSEbiv],[Nparax2 Nparabiv],Npoints);


disp(' ');
disp('The F-test for the univariate fit object GM_A_C vs. the bivariate fit object GM_A_B_C had the inputs');
disp('Nparams_univariate; Nparams_bivariate; Nsamples; SSE_univariate; SSE_bivariate of:');
disp([num2str(Nparax1) '; ' num2str(Nparabiv) '; ' num2str(Npoints) '; ' num2str(SSEx1) '; ' num2str(SSEbiv)]);
disp(' ');
disp('The F-test for the univariate fit object GM_B_C vs. the bivariate fit object GM_A_B_C had the inputs');
disp('Nparams_univariate; Nparams_bivariate; Nsamples; SSE_univariate; SSE_bivariate of:');
disp([num2str(Nparax2) '; ' num2str(Nparabiv) '; ' num2str(Npoints) '; ' num2str(SSEx2) '; ' num2str(SSEbiv)]);
disp(' ');disp(' ');disp(' ');

%% WM
fig=figure(103);
set(fig,'Position',[0 0 1000 1000])

disp(['2D WM fit for ' paramA ' and ' paramC]);
%fig bottom left
temp=zeros(10,10);
temp(1:4,8:10)=1;
subplotpos=find(temp==1);
[outputdata.WM_A_C.x,...
    outputdata.WM_A_C.y,...
    outputdata.WM_A_C.z,...
    outputdata.WM_A_C.regressionresult,...
    outputdata.WM_A_C.x_unprocessed,...
    outputdata.WM_A_C.y_unprocessed,...
    outputdata.WM_A_C.z_unprocessed,...
    outputdata.WM_A_C.GMmask,...
    outputdata.WM_A_C.WMmask] = marmplots(paramA,paramC,[],whitemask,103,subplotpos);

disp(['2D WM fit for ' paramB ' and ' paramC]);
%fig bottom right
temp=zeros(10,10);
temp(7:10,8:10)=1;
subplotpos=find(temp==1);
[outputdata.WM_B_C.x,...
    outputdata.WM_B_C.y,...
    outputdata.WM_B_C.z,...
    outputdata.WM_B_C.regressionresult,...
    outputdata.WM_B_C.x_unprocessed,...
    outputdata.WM_B_C.y_unprocessed,...
    outputdata.WM_B_C.z_unprocessed,...
    outputdata.WM_B_C.GMmask,...
    outputdata.WM_B_C.WMmask] = marmplots(paramB,paramC,[],whitemask,103,subplotpos);

disp(['3D WM fit for ' paramA ' and ' paramB ' and ' paramC]);
%fig top
temp=zeros(10,10);
temp(2:9,1:6)=1;
subplotpos=find(temp==1);
[outputdata.WM_A_B_C.x,...
    outputdata.WM_A_B_C.y,...
    outputdata.WM_A_B_C.z,...
    outputdata.WM_A_B_C.regressionresult,...
    outputdata.WM_A_B_C.x_unprocessed,...
    outputdata.WM_A_B_C.y_unprocessed,...
    outputdata.WM_A_B_C.z_unprocessed,...
    outputdata.WM_A_B_C.GMmask,...
    outputdata.WM_A_B_C.WMmask] = marmplots(paramA,paramB,paramC,whitemask,103,subplotpos);

set(findall(fig,'-property','FontSize'),'FontSize',12)

print(['Output/' paramA '_' paramB '_' paramC '_WM'],'-dpng','-r300');

Npoints = outputdata.WM_A_B_C.regressionresult.fres.NumObservations;
SSEx1 = outputdata.WM_A_C.regressionresult.fres.SSE;
Nparax1 = outputdata.WM_A_C.regressionresult.fres.NumCoeffs;
SSEx2 = outputdata.WM_B_C.regressionresult.fres.SSE;
Nparax2 = outputdata.WM_B_C.regressionresult.fres.NumCoeffs;
SSEbiv = outputdata.WM_A_B_C.regressionresult.fres.SSE;
Nparabiv = outputdata.WM_A_B_C.regressionresult.fres.NumCoeffs;

disp(' ');
outputdata.WM_A_B_C.regressionresult.fres
outputdata.WM_A_C.regressionresult.fres
outputdata.WM_B_C.regressionresult.fres
disp(' ');

disp('Testing univariate fit object WM_A_C (Fit 1) against bivariate fit object WM_A_B_C (Fit 2) with a significance threshold of 0.001:');
[Fbivx1,pbivx1]=ftest([SSEx1 SSEbiv],[Nparax1 Nparabiv],Npoints);
disp(' ');
disp('Testing univariate fit object WM_B_C (Fit 1) against bivariate fit object WM_A_B_C (Fit 2) with a significance threshold of 0.001:');
[Fbivx2,pbivx2]=ftest([SSEx2 SSEbiv],[Nparax2 Nparabiv],Npoints);


disp(' ');
disp('The F-test for the univariate fit object WM_A_C vs. the bivariate fit object WM_A_B_C had the inputs');
disp('Nparams_univariate; Nparams_bivariate; Nsamples; SSE_univariate; SSE_bivariate of:');
disp([num2str(Nparax1) '; ' num2str(Nparabiv) '; ' num2str(Npoints) '; ' num2str(SSEx1) '; ' num2str(SSEbiv)]);
disp(' ');
disp('The F-test for the univariate fit object WM_B_C vs. the bivariate fit object WM_A_B_C had the inputs');
disp('Nparams_univariate; Nparams_bivariate; Nsamples; SSE_univariate; SSE_bivariate of:');
disp([num2str(Nparax2) '; ' num2str(Nparabiv) '; ' num2str(Npoints) '; ' num2str(SSEx2) '; ' num2str(SSEbiv)]);
disp(' ');disp(' ');disp(' ');


%diary('off');